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cn . Abstract 



The article is devoted to numerical studies of atomic (metal) hydrogen with Path 



CN ! Integral Monte Carlo (PIMC) technique. The research is focused on the range of 

^ I temperatures and densities where quantum statistics effects are crucial for electrons 

and negligible for protons. In this range the equations of state are obtained as 



a dependence of internal energy and pressure on temperature and density. These 
dependences allow to detect and describe the phase transition between solid and 
hquid phases. 



td ; 1 Introduction 



One of the major recent achievements of astrophysics is the discovery of numerous exoplanetary 



^ '. systems. Almost thousand such planets have been discovered [T]. Most of them are gas giants 
^ I up to ten Jovian masses. That is the reason that attracts an increasing interest to the models 
of planetary evolution. By the current conception gas giants mainly consist of hydrogen and 
helium. So the equation of state of these elements is crucial for the models of planetary 
I> . formation and evolution. 

The detection of huge magnetic moment of the solar system gas giants has proven that they 



lO ! have liquid metal hydrogen core [2]. There are some exoplanets that are much more massive 
^ I and maybe colder than Jupiter. Because of higher pressure and less temperature their cores 
^ ■ may contain not only liquid, but also solid crystal hydrogen. The formation, evolution and 
(^ ■ properties of planets are determined by the balance of gravity and pressure, and the pressure 
in one's turn is determined by the equation of state and thermodynamical parameters of the 
planetary matter. In the cores of gas giants the prevailing substance is metal hydrogen. It can 
S^ I be described as a many-body quantum system. Its analytic analysis is extremely complicated, 
so the numerical calculations are actual in this problem. 

This article is devoted to Path Integral Monte Carlo simulation of metal hydrogen. In 
the explored range of temperatures and densities electrons form a degenerate quantum gas 
while nuclei can be examined with classical statistics, that allows to avoid fermion statistics 
problem. The parameters to be explored are internal energy and pressure and their dependence 
on temperature and density. We also focus on the phase transition between liquid and crystal 
phases. It is detected and explored in a wide range of densities. 

It should be noted that the study of metal hydrogen is important not only for astrophysics, 
but also due to the progress in diamond anvil cell experiments that have recently obtained 
crystal metal hydrogen in the laboratory [3j. 

We broadly use nuclear units in this work: kg = h = e = rrip = 1; here kg is Coulomb 
constant and rUp is proton mass. Corresponding units of principal physical quantities are: 
nuclear Bohr radius aoAr = Lj^ = 2.9 x 10"^"^ m for length, nuclear Hartree Ha = Ejq = 

1 



8.0 X 10"^^ J for energy, p^ = 3.3 x 10^^ Pa for pressure unit, pn = 7.0 x 10^^ kg/ni^ for 
density unit and T/v = 5.8 x 10^ K for temperature unit. In nuclear units electron mass is 
rUe = 5.4 X 10~^ and (electron) Bohr radius is aoe = 1/^e = 1-8 x 10^. 

We study the dependence of the atomic hydrogen properties on temperature and density, 
described by parameters /3 

P = l/ksT (1) 

and Tg (Wigner-Seitz radius) 

P = ^^ (2) 

3 s 

respectively. We simulate a finite cell of the substance, containing N^ = 128 particles. The 
properties to be evaluated in the simulation are internal energy (the sum of kinetic and potential 
energies of the particles) 

E = K + V (3) 

and pressure P. It is well known that the functions E{p,T) and P{p,T) provide a complete 
thermodynamical description of the system. In order to obtain an obvious measure of the order 
of the system we also calculate Lindemann ratio 

C^^. (4) 

Here (x^) is the particle displacement from its site in crystal lattice and i?„ is the distance to 
the nearest neighbouring particle. Lindemann ratio is used to explicitly distinguish chaotic and 
crystal phase. 

2 Model 

The Hamiltonian of atomic hydrogen is 

Hfuii = KN + K, + Vo + Ve + Vm- (5) 

Here i^jv and K^ are kinetic energies of nuclei (protons) and electrons respectively. Vq, K 
and Vint are potential energies of nuclei-nuclei, electron-electron and nuclei-electron interaction 
respectively; all three are sums of pair Coulomb interaction, for example 

«1=1 12 = 1 

There is a wide range of temperatures and densities where on the one hand electrons can be 
considered as degenerate Fermi gas and Tomas-Fermi model is applicable to them (i. e. Fermi 
statistics is of primary importance), but on the other hand protons are strongly not degenerate 
and their statistics is of no importance. On these assumptions we can deal only with protons, 
moreover we can use classical (Botzmann) statistics. The effect of taking electrons into account 
is Thomas-Fermi screening. So the effective Hamiltonian is 

Hfuii = Kn + Vn. (7) 

Here V^ is potential energy of protons with screened interaction: 

Np Ji-1 



y ^ v^ v^ exp{-ri^ij Rtf} 



n=l J2=l 



Thomas-Fermi screening length Rtf is given by 

Rtf = a—^aoeVs. (9) 

The Hamiltonian ([5]) can be reduced to ([7]) under following conditions. First, we want to 
neglect effects of nuclear forces for protons, so their separation (which is approximately r^) 
must me much greater than their size Rp. Second, we want to applicate Thomas- Fermi theory 
to electrons, that can be done if there are many electrons within screening length. This leads 
intuitively obvious restriction that nuclei separation must me less than (electron) Bohr radius. 
So, our approximation are applicable for densities corresponding 

Rp<^rs<^ aoe- (10) 

The limits in nuclear and SI units are Rp ^ 3 x 10"^ ^ 9 x 10"^^ m and aoe ~ 2 x 10^ ^ 

5 X 10~^^ m. The estimations for limiting densities ([2]) are pmin ~ 3 x 10^ kg/m^ and pmax ~ 

6 X 10^^ kg/m^. 

Next, our approximation is valid if electrons are degenerate and protons are not. The 
degeneracy temperature can be estimated as (3d ~ mr"^. So, the acceptable range of temperatures 
depends on density and it is defined as 

ruerl ^P -^r]. (11) 

The temperature limits in nuclear units are Pmin ~ 5 x 10""^ r^ and [imax ~ ^1- This leads 
following estimations at given densities (in SI units): T^in ~ 0.9p3kg~^/^m^K and T^ax ~ 
2 X lOV^kg-^/^m^K. 

3 PIMC 

3.1 Path Integral Monte Carlo 

Suppose a system, determined by coordinates x, in imaginary time. The density matrix of such 
system with Hamiltonian H at the temperature /3 is 

Pxo^xjv, = (xo|e"'^^|x^,). (12) 

Its partition function is 

^ = trp= /"cixo(xo|e-''^|xo). (13) 

Average observable A is calculated with 

{A) = ^tr(Ap) = iy"rfxo(xo|Ae-^^|xo). (14) 

To proceed to the path integral formulation, introduce the "time step" r, defined as 

1/T = I3 = NtT. (15) 

and decompose the density matrix into a product of Nt density matrices 

Pxo^xjvj Pxo^xi • • • Pxt-i^rxt ■ ■ ■ Pxjvj - 1 ->-a;jVi ) \^^) 



where each of these intermediate matrices is 

Px,_,^x. = (xi_i|e-^^|xi) = e-'\ (17) 

"Lattice action" S is defined as 

Nt 

S = J2St. (18) 

In fact the above decomposition given by Trotter formula is correct only if A^j — > oo, and for 
real simulation A'"( will be chosen large enough to eliminate the dependence of the result on it. 
Next we introduce the notation 

Nt 

I?x = J]rfxt (19) 

t=i 

and consequently the formulae ( 1T3|) and ( TT4|) can be represented as follows: 

rixe^^, (20) 



<'^) = Jv^pr = / -^j^^ (21) 

Formula (ETj) reveals the idea of Path Integral Monte Carlo. Since we have a (large enough) 
set of paths a; = xq . . . xt . . . xat^ , where the probability of the path to be included into the set 
is proportional to its "statistical weight" 

7r{x) ~ e-^("). (22) 

The average of any observable can be measured by simple (arithmetic) averaging over this set. 

3.2 Algorithms 

The way to obtain properly distributed ( 122|) paths is based on the property of Markov chains to 
converge to the limiting distribution. A sufficient condition of the convergence to the limiting 
distribution 7r(x) for the Markov chain with a transition probability V{x — )■ x') is the detailed 
balance condition: 

V{x-^x')tt{x)=V{x' ^ x)tt{x'). (23) 

The specific form oi V{x — )■ x') is not fixed, but it must be constructed carefully as it crucially 
affects the time of "thermalization" (convergence to the limiting distribution). 

A generalized Metropolis-Hastings algorithm is based on the decomposition of transition 
probability: 

n. -. .', ^ r(. -. .')Ai. -. .') + Si. - .', {l -J^yn. -. y)Ai. ^ y)] . (24, 

here 



A{x — )■ x') = min 



T{x' — )■ x)tt{x') 



(25) 



T{x — > x')tt{x) 

It satisfies the detailed balance condition for any T{x — ?■ x'). Formula ( l25l) means the following. 
First, generate a new (trial) configuration with probability T; then accept it (add it to the set) 
with probability A or reject it (return to the previous configuration and add an other copy of 



it to the set) with probabihty 1 — ^. The specific form of the algorithm is defined by the choice 
of the function T{x — )■ x'). The theoretically best choice is "heat bath": 

r{x ^ x') = 7r(x'), A{x -^ x') = 1, V{x -^ x') = 7r(x'). (26) 

Unfortunately, most probability distributions can not be generated directly fast enough, so we 
have to use a general type of the algorithm fl25|) . There are two demands to the distribution 
T{x — )■ x') : first, it must be close to 7i{x'), second, there must be an algorithm of generating 
it numerically very fast. It is rather natural to choose T(x — )■ x') as the kinetic part of the 
"statistical weight", then the acceptance probability A{x — t- x') is proportional to its potential 
part. 

Primitive algorithm is based on "sweep" when the transition from "old" configuration to 
"new" one is a try to change only one coordinate (or coordinates in the only imaginary time 
slice t). For large systems and for large number of slices it has huge autocorrelation. It means 
that "new" configurations turn out to look like "old", and it takes much time to obtain really 
statistically independent ones. This problem can be solved with the multilevel algorithm [Ij. It 
is based on fast generation of a rough approximation of the path, that increase the acceptance 
rate of the further more accurate one. 

Consider a bisection multilevel algorithm. We start from a part of the path with length 
2'^ievei slices, for example s = (xq, . . . ,X2iV;e„e!)- This part of the path is divided into levels 
Sfc. Zero level consists of the coordinates on the boundaries of the chosen part of the pass: 
So = (xo,X2iV;g„g,). They are not to be changed during the current multilevel update. The 
first level consist of the coordinates on one middle time slice Si = (xgAfj^^^j-i). The second level 
consist of two time shces S2 = {yi2^ievei-2 ,^2^ievei-i+2^ievei-'2) , etc. There are 2''~^ shces in the 
k-th level. Introduce "level action" iik^Sk) = 7rfc(so, . . . , Sk-i, Sk), which is a function of Sk and 
previous levels coordinates are parameters. Intermediate levels actions can be chosen arbitrary, 
the only requirement is that the action of the last level must be the lattice action: 

7r7v,_,(s7v,_J = 7r(s). (27) 

Then start a kind of Metropolis-Hastings algorithm with trial probability distribution 

'^k{s^.) = Tk{sQ, . . . , s,^_-^; Sk] Sk+i, . . . , Sni^^^i — )■ Sg, . . . , Si._i; s,.; Sk+i, . . . , Sni^^^i) 

and acceptance probability 



•Ak{s'k) — A.{Sq, . . . , 4_i; Sk] Sk+1, . . . , SNi^^^i -> s'q,. . . , S^_i; s'f.] Sk+l, . . . , SNi^^^i) 



mm 



1, 



It satisfies the level detailed balance condition 



(28) 



' k\^k)- — t; — T - ykKSk)- — T-, — r (29) 

7rfc_i(Sfc_i) T^k-l{Sf^_^) 

that leads to full detailed balance ( 123|) : 

T^k{Sk) = / dSk+1 . . . dsNi,,,,Tr{s). (30) 



3.3 Some Details 

Our simulation is limited in the number of particles, and consequently in the spatial size of the 
cell. We use cubic cell and periodic boundary conditions in space. The size of the cell is 



"^ttNts. (31) 

The a = x,y,z coordinate of i-th particle in the t-th time slice is denoted by x'^{t). To 
describe the configuration completely we also need "winding numbers" n"(t) = —1,0,1, that 
denote if the corresponding path "skips" from one side of the cell to another through peri- 
odic spatial boundary conditions. Potential energy of particle interaction (particles can be in 
different " copies" of the cell due to boundary conditions) is determined by their separation 



<.f "'(^) = \ J2('^m - xm + Ln")2. (32) 

\ a=l 

In the notation, described above, the lattice action corresponding to the Hamiltonian (^^ with 
potential energy ([H]) and periodic spatial boundary conditions is set as 

— Invr = S = St + Sy- (33) 

S^ - V V V ^^^^^^ ~^^^^~ ^^ + ^<^^^^' , (34) 

t=l i=l a=l 

In the case of periodic boundary conditions the trial probability density based on the kinetic 
part of the action can be represented as (skipping irrelevant indices for simplicity) 



r(x°°(t)|n(t + 1) - n{t)) ~ exp ■ ^ 

r 



x°°(t) 



x{t + 1) + x{t - 1) + L{n{t + 1) - n(t)) 



2 



2 

^ (36) 

Gaussian (it has infinite range) distribution of x°°{t) = x'{t) + Ln'{t) can be generated fast (we 
use Box-Muller transform) and allows to determine n'(t), n'{t + 1) and x'{t) due to conditions 
-L/2 < x'{t) < L/2 and n'{t + 1) - n'{t) = n{t + 1) - n{t). 

We use multilevel algorithm. Though the level action can be chosen arbitrary, there is a 
theoretically optimal choice. The action of the level should be obtained by integrating out the 
next levels coordinates in the full lattice action: 

vrfc(sfc) = / dsk+i . . . rfs7v,_iVr(s). (37) 

For our model with action fl33|) . flMl) . fl35l) it leads to a quite simple and effective algorithm. Trial 
probability distribution for each bisection is fl36|) . where the level time step is r — )■ r^ = 2^''="'='^'^r 
and winding number conserves nfc(t + l) — nk{t) = nk-i{t). This trial distribution together with 
the condition ( 1371) leads to the acceptance probability 



Akis'k) = min 



-Svisk) 



(38) 



Sv{sk) is determined by ( 135|) with the first sum only over the slices that belong to the level s^. 
T is not level time step (as it was in the kinetic part) but the real time step. 



4 Results 

The calculations were performed for following parameters. Na = 1//3 from 0.5 x 10^^ to 
4.75 X 10~^ with step 0.25 x 10~^ and for additional points 0.57 x 10"^ and 0.66 x 10^5. r^ 
was changed from 200 to 450 with step 50. These values in nuclear units correspond to SI 
values of temperature from 2.9 x 10^ K to 27.7 x 10^ K and density from 183 x 10^ kg/m^ to 
2085 X 10^ kg/m^. The lattice of calculation points will be shown in Figure [T6] (discussed later). 

4.1 Energy 

Average internal energy {E) is calculated as ([3]), taking into account (IT5l) . (lMl) . (l35l) : 



{V) = {^), (39) 

(A-) = (^ - f ). (40) 

Note that while the potential energy observable is rather intuitive, the kinetic energy one is 
quite different from intuitive (but incorrect) form. By the way in real numerical calculations 
the averaging should be done exactly as in ( HUj) . (T) = 3iVp/2r — {St) / P seems similar but 
leads to large errors because of substraction of very close large numbers. 

Figures [1] and |2] show the internal energy ii^ as a function of temperature for the densities 
183 X 10^ kg/m^ and 2085 x 10^ kg/m^ respectively. In both cases we observe a slight increase 
with increasing temperature and an acute jump at certain temperature that is associated with 
the phase transition. Figures E] and H] show the potential energy at these densities, which 
behaves similar to full energy, i. e. increases and has a jump up at the same temperatures for 
each given density. Figures |5] and [6] show the kinetic energy at these densities. Its behaviour is 
different from potential and full energy. It also increases, but jumps down at phase transition. 
At lower densities this jump vanishes and turns into a jump of the slope only. So, it looks hke 
a second-order phase transition at densities 261 x 10^ kg/m^ and lower and like a first order 
phase transition at densities 618 x 10^ kg/m^ and higher. 

We can see that the properties of the system depend on density much stronger than on 
temperature. Correspondingly, the internal energy almost totally consist of potential energy 
determined by the distance between protons i. e. by density. In spite of this fact, the jumps 
of both parts of energy at the phase transition are of close magnitudes. In order to extract the 
main term we introduce Vqk - potential energy of "ideal zero temperature" crystal. It means 
that the particles in this crystal are exactly in the sites of its bcc (body-centric cubic) lattice 
(in all time slices). Vqk depends only on density and this dependence is shown in the Figure 
[71 We substract this zero energy from full and potential energy in order to extract non-trivial 
terms. It turns out that substracted full and potential energy and kinetic energy are of the 
same magnitude; their dependences on density and temperature also have close magnitudes. 
The substracted full internal, substracted potential and kinetic energies for a range of densities 
between 183 x 10'^ kg/m^ and 2085 x 10^ kg/m^ are shown in figures [HI [9l and [TOl respectively. 

4.2 Pressure. Equation of state 

The observable for pressure is 

(^) = ^|W-^(>:^r.,)|. (41) 
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Figure 1: E{T) at p = 183 x 10^ kg/m^ (r, = 450) 
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Figure 2: E{T) at p = 2085 x 10^ kg/m^ (r^ = 200) 

Figures [11] and [12] show the temperature dependence of pressure at above mentioned den- 
sities. It has a jump at the same temperatures as energy that proves the existence of phase 
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Figure 3: V{T) at p = 183 x 10^ kg/m^ (r, = 450) 
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Figure 4: V{T) at p = 2085 x 10^ kg/m=^ (r^ = 200) 



transition. The expression ( 14T|) allows to determine Pok similar to Vqk and perform similar sub- 
straction procedure. Figure (fT3l) shows the dependence of Pqk on density and Figure [H] shows 
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Figure 5: K{T) at p = 183 x 10^ kg/m^ (r, = 450) 
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Figure 6: K{T) at p = 2085 x 10^ kg/m^ (r, = 200) 

the substracted pressure for all the range of explored densities. Similar to energy, pressure 
mainly depends on density and quite slightly changes with temperature. In fact it is just what 
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Figure 8: E{T) — Vqk at different densities 

should be expected for condensed matter. Tlie disadvantage of tliis property is a trouble with 
thermodynamical calculations due to orders of magnitude difference between partial derivatives 



11 



15- 



o 

^ 10- 



> 

I 

> 



5- 




"■ — ' — r 
5 



"■ — ' — r~ 

10 



— I — 1 — 1 — 1 — 1 — |— 

15 20 



"■ — ' — r~ 
25 



T, x10 K 



Figure 9: V(T) — Vqk at different densities 
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Figure 10: K{T) at different densities 

by density and by temperature. An intuitive illustration can be seen in Figure [Ml Formally 
the function P{p,T) allows to determine isobars, but the resolution of experimental data is 
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insufficient despite of quite large number of points. We can only say that isobars are some lines 
close to lines of constant density, but having some little unknown slope. This problems can be 
solved in different ways, but they are not to be discussed here. 
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Figure 11: P(r) at p = 183 x 10^ kg/m^ 
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4.3 Phase transition 

The Lindemann ratio (j4]) is a good measure of disorder of the lattice, so it is extremely useful 
and obvious to detect the phase transition, where the order totally vanishes. Figure [15] shows 
Lindemann ratio for all the range of explored densities and temperatures. The phase transition 
is clearly seen here. Note that while the plateau in solid phase (left bottom) gives some 
physical information, the plateau in liquid phase (right top) is due to finite volume effects and 
it is determined only by the volume. 

The position of the phase transition is determined quite accurate, so we can draw the phase 
plane for metal hydrogen. It it shown in the Figure [121 

Here we have to describe some important details. As we know, the PIMC observables are 
averages over a set of thermalized path. To get this set we start with any path and perform a 
Markov chain procedure called thermalization. Sometime we start to get thermal equilibrium 
paths, but we do not know how soon it will be. It is well known that models of systems near 
a phase transition are usually difficult to be thermalized over the transition. For example we 
start our simulation with ideal "zero temperature" crystal lattice (solid state). During the 
calculations it thermalizes quite fast to some other solid state that seems stable. It takes quite 
much calculation time to receive true physical paths. An example is shown in Figure [T71 This 
is the main obstacle to determine the position of the phase transition more accurate. Moreover, 
it turns out that thermalization from liquid to solid state takes so much time that it hardly 
ever can be performed in moderate time. So, the position of the phase transition is formally the 
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Figure 12: P{T) at p = 2085 x 10^ kg/m^ (r, = 200) 
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Figure 13: P{T) at different densities 

upper limit. The lower limit must formally be determined with a series of simulations starting 
from "liquid" path. But it is not expected to differ much from the upper limit that we received. 
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Figure 14: P(T) — Pqk at different densities 
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Figure 15: Lindemann ratio at different densities 



15 




p, x10^ kg/m^ 



Figure 16: Phase plane 



-liquid 
- crystal 



c 
E 

T3 



1,0- 



0,8- 



0,6- 



0,4- 



0,2- 



0,0- 




20 



40 



T" 



"T 



"T 



"T 



60 



80 100 120 140 



Figure 17: Lindemann ratio thermalization at p = 2085 x 10^ kg/m^, T = 13.1 kK (blue) and 
T = 14.5 kK (red) 
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5 Conclusions 

Path integral Monte Carlo technique was implemented to simulate atomic metal hydrogen from 
the first principles. Its thermodynamical properties were explored in a wide area of the phase 
plane. Numerical equations of state were obtained. The phase transition between liquid and 
solid crystal phases was detected and explored. 

The principal thermodynamic parameters: temperature, density, pressure and energy were 
set or measured, but entropy was not. That will be the object of our following studies. The 
algorithms of obtaining entropy and adiabats are a little bit more complicated than for isoterms 
for example, because entropy can not be measured as a PIMC observable. So we have to solve 
differential equations derived from thermodynamics. Formally they give all information about 
the system since we know E{p, T) and P{p, T), but it is not trivial to get the numerical results. 
As it was mentioned, the lattice of calculation points in temperature and extremely in density 
must include close points in a large range that means much calculations. On the one hand we 
want to explore a wide range. On the other hand, the points must be close enough to allow the 
calculation of derivatives. We also plan to develop an alternative way of derivatives calculation 
based on constructing observables for them. 

An other problem to be explored is to perform the thermalization from "liquid" to crystal 
solid state in order to determine the lower limit for the phase transition. We expect that it can 
be done much faster starting from two phase system. 
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